# Markdown
        library(bookdown)
    # General data analysis and transformation
        library(readr)
        library(readxl)
        library(dplyr)
        library(janitor)
        # library(gtools)
        library(stringr)
        library(units)
        library(tidyr)
        library(rlang)
    # API-related
        library(jsonlite)
        library(urltools)
    # Mapping and GIS operations
        library(sf)
        library(leaflet)
        library(htmlwidgets)
        library(geojsonsf)
        library(rmapshaper)
        library(tmap)
        library(tmaptools)
        library(ceramic)
    # workflow
        library(here)
    # plotting
        library(ggplot2)
        library(forcats)
# Get list of CES 3 variable names and plain text names
    ces3_variables <- readr::read_csv(here('data_prepared', 'ces_names.csv'))

ces_measures <- ces3_variables %>% 
    filter(group == 'pollution burden', 
           subgroup == 'exposures', 
           type == 'percentile')

Introduction

This document attempts to do some quantitative geospatial analysis to look at correlations between redline mapping and indicators of water quality, public health, and facility information. It computes the area weighed average percentile for different CalEnviroScreen 3.0 (CES) indicators by HOLC rating for a given city and statewide. This document considers indicators in the ‘Exposures’ category of the CES Pollution Burden indicators (factsheet, indicators page), including:

  • Ozone Percentile
  • PM 2.5 Percentile
  • Diesel PM Percentile
  • Drinking Water Percentile
  • Pesticides Percentile
  • Tox Releases Percentile
  • Traffic Percentile

The CES polygons and the Redline polygons have different coverages, as shown in the maps below. In the first map for each indicator, the left pane of each map shows the HOLC rated areas, the center pane shows the CES indicator scores (at the census tract level) with the Redlined areas superimposed on top, and the right pane shows the overlapping portions of the CES polygons and the redlined areas (border colors represent HOLC rating).

    # load data
        # Redline Polygons
            redline_polygons <- read_rds(here('data_prepared', 'redline_polygons.RDS'))
                # st_crs(redline_polygons) # to check the reference system
            # modify the 'city' column to better keep track when joining with other datasets
                redline_polygons <- redline_polygons %>% rename('redline_city' = 'city')
            # add the area of each redline polygon
                redline_polygons <- redline_polygons %>% mutate(redline_poly_area = st_area(.))
    
        # CES 3 Polygons
            ces_3_poly <- read_sf(here('data_raw', 'CalEnviroScreen3', 'CESJune2018Update_SHP','CES3June2018Update.shp'))
                # st_crs(ces_3_poly) # to check the reference system
            # modify the 'City' column to better keep track when joining with other datasets
                ces_3_poly <- ces_3_poly %>% rename('CES_City' = 'City')
                
        # Transform both to UTM Zone 10N
            redline_polygons <- redline_polygons %>% st_transform(crs = 26910)
            ces_3_poly <- ces_3_poly %>% st_transform(crs = 26910)
            
    # create a map to check that the data loaded correctly
        # tm_shape(ces_3_poly) + tm_borders() + tm_shape(redline_polygons) + tm_borders(col = 'red')
    # get a basemap (See: https://github.com/mtennekes/tmap/issues/185)
        Sys.setenv(MAPBOX_API_KEY = "pk.eyJ1IjoiZGFsdGFyZSIsImEiOiJjazZqcm5ocHgwMG84M2tteGc5eWd4YXNmIn0.nHiYenoq-wFWwH3V4vtPXA")
        background <- cc_location(ces_3_poly %>% filter(CES_City == 'Sacramento'), verbose = FALSE)
    # create map
        map_redline <- tm_shape(background) + 
            tm_rgb(alpha = 0.5) + 
            tm_shape(redline_polygons %>% filter(redline_city == 'Sacramento'), is.master = TRUE) + 
            tm_polygons('holc_grade', palette = c('green', 'blue', 'yellow', 'red'), title = 'HOLC Grade') + 
            tm_layout(legend.bg.color = 'white', legend.bg.alpha = 0.7, legend.frame = TRUE)
    # get the intersection of the CES polygons and redline polygons
        ces3_redline_intersection <- st_intersection(ces_3_poly, redline_polygons)
    # inspect
        # glimpse(ces3_redline_intersection)

# add the area of each clipped polygon to the data frame
        ces3_redline_intersection <- ces3_redline_intersection %>% 
          mutate(clipped_area = st_area(.))
            # glimpse(ces3_redline_intersection)
    ces_poly_map <- function(ces_id, ces_measure_name) {
        # plot the CES polygons
        map_ces <- tm_shape(background) + 
            tm_rgb(alpha = 0.5) + 
            tm_shape(ces_3_poly %>% filter(CES_City == 'Sacramento' | CES_City == 'West Sacramento')) + 
            tm_polygons(col = ces_id, title = ces_measure_name, border.alpha = 1) + 
            tm_shape(ces3_redline_intersection %>% filter(redline_city == 'Sacramento'), is.master = TRUE) +
            tm_borders(lwd = 0, alpha = 0) +
            # tm_shape(redline_polygons %>% filter(redline_city == 'Sacramento')) +
            # tm_borders(lwd = 2, col = 'black') +
            tm_shape(redline_polygons %>% 
                       filter(redline_city == 'Sacramento', holc_grade == 'D')) + 
            tm_borders(lwd = 2, col = 'red') + 
            tm_shape(redline_polygons %>% 
                       filter(redline_city == 'Sacramento', holc_grade == 'C')) + 
            tm_borders(lwd = 2, col = 'gold2') + 
            tm_shape(redline_polygons %>% 
                       filter(redline_city == 'Sacramento', holc_grade == 'B')) + 
            tm_borders(lwd = 2, col = 'blue') + 
            tm_shape(redline_polygons %>% 
                       filter(redline_city == 'Sacramento', holc_grade == 'A')) + 
            tm_borders(lwd = 2, col = 'green') +
            # tm_text('holc_grade', size = 0.5) + 
            tm_layout(legend.bg.color = 'white', legend.bg.alpha = 0.7, legend.frame = TRUE)
        map_ces
        
    }
    ces_redline_overlap_map <- function(ces_id, ces_measure_name) {
        map_overlap <- tm_shape(background) + 
            tm_rgb(alpha = 0.5) + 
            tm_shape(ces_3_poly %>% filter(CES_City == 'Sacramento')) + 
            tm_borders(lwd = 1) + 
            tm_shape(ces3_redline_intersection %>% filter(redline_city == 'Sacramento'), is.master = TRUE) + 
            tm_polygons(col = ces_id, title = ces_measure_name, border.alpha = 0) + 
            # tm_shape(redline_polygons %>% filter(redline_city == 'Sacramento')) + 
            # tm_borders(lwd = 2, col = 'black') + 
            tm_shape(redline_polygons %>% 
                       filter(redline_city == 'Sacramento', holc_grade == 'D')) + 
            tm_borders(lwd = 2, col = 'red') + 
            tm_shape(redline_polygons %>% 
                       filter(redline_city == 'Sacramento', holc_grade == 'C')) + 
            tm_borders(lwd = 2, col = 'gold2') + 
            tm_shape(redline_polygons %>% 
                       filter(redline_city == 'Sacramento', holc_grade == 'B')) + 
            tm_borders(lwd = 2, col = 'blue') + 
            tm_shape(redline_polygons %>% 
                       filter(redline_city == 'Sacramento', holc_grade == 'A')) + 
            tm_borders(lwd = 2, col = 'green') + 
            #tm_text('holc_grade', size = 0.5) + 
            tm_layout(legend.bg.color = 'white', legend.bg.alpha = 0.7, legend.frame = TRUE)
        map_overlap
    }
    # FACET MAP
    ces_redline_facet_map <- function(ces_id, ces_measure_name) {
        map_facets <- tm_shape(ces_3_poly %>% filter(CES_City == 'Sacramento')) +
            tm_borders(lwd = 1) +
            tm_shape(background) +
            tm_rgb(alpha = 0.5) +
            tm_shape(ces3_redline_intersection %>% 
                         filter(redline_city == 'Sacramento'), is.master = TRUE) +
            # tm_shape(ces3_redline_intersection %>% filter(redline_city == 'Sacramento')) +
            tm_facets(by = 'holc_grade', free.coords = FALSE) +
            tm_polygons(col = ces_id, title = ces_measure_name, border.alpha = 0) +
            tm_shape(redline_polygons %>% filter(redline_city == 'Sacramento')) + 
            tm_facets(by = 'holc_grade') +
            tm_borders(lwd = 2, col = 'black')  +
            # # D
            # tm_shape(redline_polygons %>%
            #            filter(redline_city == 'Sacramento', holc_grade == 'D')) +
            # tm_facets(by = 'holc_grade') +
            # tm_borders(lwd = 2, col = 'red') +
            # # C
            # tm_shape(redline_polygons %>% 
            #            filter(redline_city == 'Sacramento', holc_grade == 'C')) + 
            # tm_borders(lwd = 2, col = 'gold2') + 
            # tm_facets(by = 'holc_grade') +
            # # B
            # tm_shape(redline_polygons %>% 
            #            filter(redline_city == 'Sacramento', holc_grade == 'B')) + 
            # tm_borders(lwd = 2, col = 'blue') + 
            # tm_facets(by = 'holc_grade') +
            # # A
            # tm_shape(redline_polygons %>% 
            #            filter(redline_city == 'Sacramento', holc_grade == 'A')) + 
            # tm_borders(lwd = 2, col = 'green') + 
            # tm_facets(by = 'holc_grade') +
            #
            tm_legend(bg.color= 'white', frame = TRUE)

        map_facets
    }
ces_scores_plot <- function(ces_id, ces_measure_name) {
    
    # group by city and holc rating, then compute the weighted average score for each
        ces3_redline_grouped <- ces3_redline_intersection %>% 
            st_drop_geometry() %>%
            mutate(area_x_score = clipped_area * (!!as.name(ces_id))) %>% 
            group_by(redline_city, holc_grade) %>% 
            summarize(total_area = sum(clipped_area), 
                      area_x_score_total = sum(area_x_score)) %>% 
            mutate(weighted_score = area_x_score_total / total_area) %>% 
            drop_units()

    ces3_redline_grouped_statewide <- ces3_redline_grouped %>%
            group_by(holc_grade) %>%
            summarize(total_area = sum(total_area),
                      area_x_score_total = sum(area_x_score_total)) %>%
            mutate(weighted_score = area_x_score_total / total_area) %>%
            mutate(redline_city = '**Statewide**') %>%
            drop_units()

    ces3_redline_grouped_city_state <- bind_rows(ces3_redline_grouped, ces3_redline_grouped_statewide)

    # plot the results by city and holc rating, and include statewide summary
        g_city_state <- ggplot(ces3_redline_grouped_city_state) +
            aes(x = fct_relevel(fct_reorder(redline_city, weighted_score), 
                                c('**Statewide**'), 
                                after = 0), 
                y = weighted_score, 
                fill = fct_rev(holc_grade)) +
            geom_bar(stat = 'identity', position = 'dodge') +
            geom_vline(xintercept = 1.52, size = 0.5, color = 'grey50') +
            scale_fill_manual(values = rev(c('green', 'blue', 'yellow', 'red'))) +
            guides(fill = guide_legend(reverse = TRUE)) +
            labs(fill = 'HOLC Rating', 
                 title = ces_measure_name,
                 x = NULL, 
                 y = "Area Weighted Average Score By City") +
            coord_flip()

        g_city_state
}

Pollution Burden - Exposures

Ozone Percentile

    i_map <- 1
    ces_measure_name <- ces_measures$name[i_map]
    ces_id <- ces_measures$id[i_map]
    
    map_ces <- ces_poly_map(ces_id, ces_measure_name)
    map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    
    tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)

    ces_redline_facet_map(ces_id, ces_measure_name)

    ces_scores_plot(ces_id, ces_measure_name)

PM 2.5 Percentile

    i_map <- 2
    ces_measure_name <- ces_measures$name[i_map]
    ces_id <- ces_measures$id[i_map]
    
    map_ces <- ces_poly_map(ces_id, ces_measure_name)
    map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    
    tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)

    ces_redline_facet_map(ces_id, ces_measure_name)

    ces_scores_plot(ces_id, ces_measure_name)

Diesel PM Percentile

    i_map <- 3
    ces_measure_name <- ces_measures$name[i_map]
    ces_id <- ces_measures$id[i_map]
    
    map_ces <- ces_poly_map(ces_id, ces_measure_name)
    map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    
    tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)

    ces_redline_facet_map(ces_id, ces_measure_name)

    ces_scores_plot(ces_id, ces_measure_name)

Drinking Water Percentile

    i_map <- 4
    ces_measure_name <- ces_measures$name[i_map]
    ces_id <- ces_measures$id[i_map]
    
    map_ces <- ces_poly_map(ces_id, ces_measure_name)
    map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    
    tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)

    ces_redline_facet_map(ces_id, ces_measure_name)

    ces_scores_plot(ces_id, ces_measure_name)

Pesticides Percentile

    i_map <- 5
    ces_measure_name <- ces_measures$name[i_map]
    ces_id <- ces_measures$id[i_map]
    
    map_ces <- ces_poly_map(ces_id, ces_measure_name)
    map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    
    tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)

    ces_redline_facet_map(ces_id, ces_measure_name)

    ces_scores_plot(ces_id, ces_measure_name)

Tox Releases Percentile

    i_map <- 6
    ces_measure_name <- ces_measures$name[i_map]
    ces_id <- ces_measures$id[i_map]
    
    map_ces <- ces_poly_map(ces_id, ces_measure_name)
    map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    
    tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)

    ces_redline_facet_map(ces_id, ces_measure_name)

    ces_scores_plot(ces_id, ces_measure_name)

Traffic Percentile

    i_map <- 7
    ces_measure_name <- ces_measures$name[i_map]
    ces_id <- ces_measures$id[i_map]
    
    map_ces <- ces_poly_map(ces_id, ces_measure_name)
    map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
    
    tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)

    ces_redline_facet_map(ces_id, ces_measure_name)

    ces_scores_plot(ces_id, ces_measure_name)